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1. Introduction 

During the last decade, independent observational data such as type-la Supernovae (SNIa) 
m, ^ and Cosmic Microwave Background (CMB) (Bj, ^ and Baryonic Acoustic Oscilla- 
tions It], ^, suggest that two dark components govern the dynamics of the Universe. They 
are the Dark Matter (DM), thought to be the main responsible for structure formation, 
and a non-zero cosmological constant A (see, e.g. ref. ||l^) or a dynamical dark energy 
(DE) component, that is supposed to drive the observed cosmic acceleration |11, 12|. 



The standard cosmological model ACDM provides a good fit to observations, but it 
describes the Universe by means of two unknown components, which represent 96% of the 
total energy density. However, it should be recognised that, while some form of Cold Dark 
Matter (CDM) is independently expected to exist within any modification of the Standard 
Model of high energy physics, the really compelling reason to postulate DE has been the 
discovery that the Universe is experiencing a phase of accelerated expansion. This could 
be theoretically unsatisfactory and, recently, alternative models have been proposed with 
respect to the class of DE models. In particular, many of them are based on modifications 
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of gravity at large distances. For example, scalar-tensor theories |TJ, |l5|, 16, 17, 18, 19 1, 
Brane- World models (see e.g the review [^] and refs. therein), Galileon models [pl| , 22, 
23, U, 25, 26, 27], Gauss-Bonnet gravity |28, ^] and other scenarios (see, for example, 
the review |11] and refs. therein). 

An interesting class of modified gravity models is represented by f{R) theories (see 
the reviews 32, 



p4\ , pq , pq , p7| ] and refs. therein ), whose Lagrangian density is 
simply defined by an arbitrary function of the Ricci scalar R. These Lagrangians were 
proposed for the first time in connection with Infiation in the early Universe |61] and, only 
recently, they have been used in the context of DE models, to explain the present-day 
cosmic acceleration (see e.g. [ 



63, B4 



In general, there are two ways to study these models |65]: 1) specifying directly the 
type of Lagrangian that satisfies cosmological and local gravity constraints (e.g. []66| , |67| 
and see also [^]) 2) through a parametrized post-Newtonian framework or similar (see e.g. 
Refs. [|^, 69, |7l|, ^ |7^, 74, 75|), to provide a scale-dependent parameterization of 
cosmological perturbations (see also [7^]). 

Recently, it was found that, in viable f{R) models (see e.g. |77, |6^, 



79, 80, 81 



82, 33, 56, pj]) when the super horizon long- wavelength limit is taken (i.e. when spatial 



gradients may be neglected in the equations), the evolution of metric and density perturba- 
tions of a Robertson- Walker background can be described by the Friedmann equation and 
energy-momentum conservation p^, [ 



j, 74 1 . In addition, the viable f{R) models 

usually need to be close to the ACDM model during the matter dominated epoch in order 
to satisfy cosmological constraints, local gravity constraints (see 



34, 84, 86, 



and 



galactic constraints [S7|. Motivated by these results, in this paper, we propose an ana- 
lytical method which allows to describe both the background evolution and cosmological 
perturbations in f{R) modified gravity models. This approach is completely different from 
previous ones. Indeed, through this analytic technique we can determine quantitatively the 
deviation in the behavior of cosmological perturbations between a given f{R) model and a 
ACDM reference model. Moreover, our treatment is general in that all the results depend 
only on the initial conditions that characterize the type of f{R) model we are studying. 
Finally, our approach allows to study structure formation in these models from the largest 
scales, of the order of the Hubble horizon, down to scales deeply inside the Hubble radius, 
without employing the so-called "quasi-static" approximation (see also ||8^, |8^] ) . Although 
we restrict our analysis here to linear perturbations, the technique can be applied to lin- 
ear and weakly non-linear scales and can be extended to any perturbative order (see also 
|9C, 91|). Of course, in order to study these models at non linear scales another approach 
must be adopted (for example, see |34]). 

The rest of the paper is organised as follows. In Section ^ we present the basic equations 
describing the background and the perturbative evolution of a generic f{R) gravity model. 
In Section ^ we describe our approach, analyzing both the background evolution and the 
first-order perturbation equations in the synchronous gauge and, in Section ^ we introduce 
the evolution of cosmological perturbations in the Poisson gauge. In Section ^ we recover 

'See also ||, |ll g, 0, ||, |l|, g, g g, g, g for other 

recent contributions to this area. 
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the "quasi-static" approximation and Section |6| we compute the matter power spectrum 
with our approach. Section ^ is devoted to our conclusions. Appendix |^ and Appendix ^ 
are devoted to the evolution of linear perturbations in the reference ACDM model, in the 
synchronous and Poisson gauge, respectively. 

Throughout the paper we use G = c = 1 units and the (— , +, +, +) signature for the 
metric. Greek indices run over {0, 1, 2, 3}, denoting space-time coordinates, whereas Latin 
indices run over {1,2,3}, labelling spatial coordinates. 



2. Preliminary relations 

Let us consider the following action in f{R) gravity (see e. g. |61, 62, B^, 64]; see also the 



reviews [31, 34, 35 1 and refs. therein): 

S = + = / d'xV^fiR) + I d'xV^C^[^^,g^.] . (2.1) 

In this case f{R) is a general function of the Ricci scalar, R, and £m is the matter La- 
grangian and are the matter fields. Defining cp = df /dR, S^'^^^ can be cast in the form 



of Brans-Dicke (DB) theory |13] with a potential for the scalar field (p 1 31, 33, 35]. 
In particular. 



^(GR) 



1 

16^ 



d^x^ m - y(</.)] , 



(2.2) 



where V{(j)) = R(p — fiR)- The field equations obtained from varying the action with 
respect to g^u are 



where Gfj,^ = R^^ - {l/2)Rg^ 



1, 



SttT/, = V^V.cp - U(pg^, - -V{(t>)g^, 
and is the energy-momentum tensor of a perfect fluid, i.e. 

2 



-g ^9' 



(Pm + Pm)Uf,U„ + p 



(2.3) 



(2.4) 



(2.5) 



The vector is the fluid rest-frame four- velocity, p is the energy density and p the isotropic 
pressure. In this case the equation of motion of the scalar field is simply 



(2.6) 



Finally by taking the trace of Eq. ( p.3p and using Eq. (2.6), we obtain the dynamics of the 
scalar field for a given matter source 



□</> + -( 2V 



.dV 



Svr, 



(2.7) 



where T"^ = g'^'T^. 
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When considering the background cosmological evolution, we take the metric to be of 
the flat Friedmann-Lemaitre-Robertson- Walker (FLRW) form, ds^ = a?'{'q){—drf + dx^), 
where r] is the comoving time and a(r/) the scale factor. Then the gravitational field 
equations, the scalar field equation of motion and the continuity equation become 

AH" = ^pL°) + {a- ^)n' + {A- ^yn + ^v{y.) , (2.8) 

A {2%' + H") = -Svra^pW + {A - ^) {2%' + U^) + {A - ^)" + {A - ^)'U + y , 

(2.9) 

^' + ^^ = ^^, (2.10) 
pL°^' + 3^(pL°^+pS^) = 0, (2.11) 

where ^ is a suitable constant (for its physical meaning see Ref. p^), ^ = a! ja and 
primes indicate differentiation w.r.t. r\. Note that (\)^^^ = ip, and^m^ are respectively 



the background scalar field, matter pressure and matter energy density. From Section 3.1 
onwards we will set A = 1. 

Now taking into account only the non-relativistic matter, i.e. = 0, let us consider 
the metric of a flat FLRW Universe with small perturbations. In particular we want to 
study the evolution of the scalar perturbations in the synchronous gauge (see also Ref. 
|[78|] ). In this case the line-element is written in the form ds^ = a?'{—drf + [(1 — 2il)'^^^)5ij + 
DijX^^^)\dx'^ dx^} ^ where Dij = didj — (l/3)(5ijV^ is a trace-free operator. Let us allow for 
small inhomogeneities of the scalar field, (p{r],x.) = ip{r]) + (p^-^\r] , x) . Perturbing Eq. (p.3|), 
we get 



i)A _ JO) / 




5« 

3 - y v? ■ \^ v? J 

+3«^_3«^^^-X^ + 5!(^-^i:)^, (2.12, 

if V V 2 \ (/3 Alp J tp 

2V,«' + Vx«' = ^-?^^, (2.13) 

3 ip ip 

"+(2n + ^) x(^) ' + 2^(1) + ^V^x^^) = 2^ , (2.14) 
\ ipj 3 (p 

where 6^^^ = (pm^ — Pm^)/p[n ■ From Eq. ( p.7D , we obtain 

(2.15) 
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Finally, at the linear order, Eqs. ( |2.6D and ( 2.11| ) become 



-6V(^) " - 18?^V^'^ ' + 4V2^(^) + ^V^V^x^i) = ci^^^^''^ , (2.16) 
5(i)' = 3^W'. (2.17) 
3. Description of the analytical approach 

Motivated by the fact that, in viable f{R) models |3^], 1) the cosmic scale- factor 0(77) can 
be described by the Friedmann equation with background expansion of the Universe close 
to the ACDM model and 2) gravity can be described by a classical four-dimensional metric 



theory having a well-defined infrared limit |7^, 68, 78, pq, 79, p5[ (i.e. in the modified 



gravity models, the long-wavelength perturbations obey the same constraints as they do in 



general relativity ||68|, |69|, 66, |6^]), let us make the following ansatz: 



i) at the background level, the scalar field can be described in the following way 

ip = ip + eA^ + ... , (3.1) 

where we introduced a suitable perturbative parameter |e| <^ 1; 

ii) at the first order level 

^(1) = ^(1.0) + ^ (3.2) 

^(1) =^(1,0) +,^^(1,1) + . _ (3.3) 

5(1) = 5(1,0) +,^^(1,1) + ...^ (3.4) 



and 



^2 



^(1) = ^(1,0) + ,^^(1,1) + 1.^^(1,2) + (3_5) 



where in the double superscript (i, j) the left index i refers to the order in standard pertur- 
bation theory (e.g. i = 1 denotes linear theory), while the right index j refers to the order 
of the deviation of our modified gravity model w.r.t. ACDM. In particular, setting^ 
(p = 1, A = 1 and ^(^'O) = we recover the ACDM model when e — >• 0. In other words, 
all the equations previously obtained will be also iteratively perturbed by the parameter 
e. This approach allows us to understand the relevance of the various terms that come 
from modified gravity in the description of the large scale structure of the Universe. Let us 
stress that this prescription is completely general and can be used to describe any viable 
f{R) theory on cosmological scales. 

In the next two subsections we will analyze in detail with this technique the background 
and the cosmological perturbations. Specifically we will expand all the equations in section 
|2| up to order e. 



^If we consider (p non-constant then the background could change completely. Indeed, as we will see 



from Eqs. (3.6) and (3.8), when (fi = 1 and A — 1, the value of "H does not depend on if and, obviously, on 



e^. Therefore the background is the same of the ACDM model for ^ 1. 
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3.1 FLRW background equations 

In this case, by construction, V{(f) = Vq is a constant (specifically, it is of the order of the 
cosmological constant), if = I, A = I and is unperturbed with respect to e. In this 
case from Eq. ( |2.8| ) we get 

= ^/^'m + y Vb , (3.6) 
ne + n^^ = ^V,, (3.7) 

where we have defined the expansion of V{ip) with respect to e as V{ip) = Vq + £^1(^(77)). 
Instead, from Eq. ( p.9D we obtain 

{2'H' + n^) = ^Vo, (3.8) 

i" + m' +{2'H' + 'H'^)i = —Vi. (3.9) 
We note immediately that we retrieve the field equations of general relativity at the back- 



ground level. Now from Eqs. (3/7) and ( ^.gj ) we find the equation of motion for ^ 

- 2%^' + 2{n' -n^)i = Q . (3.10) 

Notice that this field equation describes the dynamics of ^ as a free field, without any 
source term proportional to the matter fields. 

At this point, knowing that the Cauchy problem is well-formulated for metric f{R) 
gravity in the presence of matter , the initial conditions for ^ and S^' are crucial because 
they characterize the type of f{R) model that we are studying and the validity of the 
approach that we are using^. Specifically, assuming that for rj < r^recj where r^rec is some 
epoch when the Universe is matter dominated and radiation is negligible (usually at the 
recombination epoch), f{R) gravity should be the same as general relativity. Thus we can 
choose at recombination ^(??rec) = and |'^(??rec)'| ^ 1- In other words, we have to set 
both e and |C(^rec)'| "small enough" to allow only at late times a different solution of the 
perturbation evolution w.r.t. a ACDM model. 

Let us make another comment. For the stability at high curvature of perturbations of 



the scalar field the sign of e and ^' is crucial. Indeed, from Eq. ( 2.15 ) we must assume that 



(aV3) [(p{d'^V/dip^) - (dV/dtp)] > |7^, H, |7|, || . Now, taking into account that 
j2t^ f;/ 2 (0)/ (0) 

(3.11) 



for e — )• 0, we have to impose that e(ad^/da) < 0. Moreover, analyzing Eq. (3.9) we note 
that if ^' < 0(> 0) then ,^ < (> 0). Therefore, assuming for simplicity e > we have to 
set (ad^/da) < 0. 



This translates to imposing the conditions on df/dR and d^f/dR^ (or B defined by at a specific 
time (for example, at recombination or today). In particular, let us look at Eq. (3) of Ref. (see also 
Js?!]). In this paper the authors need simply to set two parameters, namely fug and n in order to study one 
class of viable f{R) models (see Section 
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Now, in order to better understand the behaviour of ^, let us rewrite Eq. ( 3.10 ) in 
a shghtly different way. Indeed, using Eqs. ( |3.6[ ) and ( |3.8| ), defining Vq = lOvrpA and 



(0) / 



(Qom/^OA), where pS^ = p|,°^c 
d^e 3 u/a 



and ao 



1, we obtain 



da^ 2 + da + 



e = o. 



(3.12) 




Figure 1: Illustrative plot of — ^, as a function of a, with (fiom/^^oA) = 3/7. The lines, from 
short to long dashes, correspond to (d^/da) (orec) — — fO^^, —10^^ respectively; the black solid line 
corresponds to (d^/da) (orcc) = — 10~^. 



In Figs. 1^ and |2| we show the evolution of {—^) and (— d^/da) as a function of the 
scale-factor a. In particular we note that the value of |^| grows rapidly after recombination 
and if the initial conditions are not sufficiently small, in future this effect could lead to 
an "explosive phenomenon" for these models |^5|. It is obvious that in that regime our 
perturbative technique cannot be used. 

Finally, it is useful to compare our formalism to the existing literature by considering 
the dimensionless quantity that quantifies the deviation from the ACDM reference model 
introduced in Ref. (see also ||0|, |3|]) 



B 



R' 



n 



m- 



where [||, g] 



In our case it becomes 



m 



B 



R 



R {%' -W^) ' 
{d'^f/dR^) {dV/d'f) 



(df/dR) ip{(PV/dip'^) 



dln(l + eg) 
da 



(3.13) 



(3.14) 



(3.15) 



*ao the value of the scalar factor today. 
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-d^/da 




0.0 



0.2 



0.4 0.6 

a 



0.8 



1.0 



Figure 2: Illustrative plot of (— d^/da), as a function of a, with (fiom/i^o a) — 3/7. The lines, from 
short to long dashes, correspond to (d^/da) (flrec) = — 10~^, — 10~^ respectively; the black solid line 
corresponds to (d^/da) (orcc) = —10^''. 



We see immediately that, imposing < and e > 0, we get the corresponding stability 
condition B > found in [77|. 



3.2 First-order perturbation equations 

In this section we analyse in detail the field perturbation equations through this iterative 



technique. From Eq. ( 2.12 ) we get 

6?^V^''°' ' - (2VV(''°^ + ^V2v2;t(i'0)) + S^VS^'^^'-"^ = , (3.16) 



(3.17) 



From Eq. (|1|) we obtain 



3 



(3.18) 
(3.19) 
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while from Eq. ( ^.14 ) 



.(1,0) // + 2-Hx(''°) ' + 2^^^'°) + - V^x^^'O) = , 



(3.20) 



^a^) " + 2nx^^''^ ' + 2V'(i'i) + - V^x^'''^ = 2e(''') - ' ■ (3.21) 

3 



Moreover, from Eq. ( 2.15 ), we get 



STra^ 



(3.22) 
(3.23) 



while, from Eq. ( p. 16 ), 



6^.(1-0) " - IST^V'^''") ' + 4V2V'^''°^ + -V^V^x^i'O) = I - 2?^3) ^(i,i) ^ (3 24) 

- eV'^^'i) " - 18?^V^'''^ ' + 4VV^'''^ + - V^V^x^^'^) = 4 (•W" - 2^') 4^''') • (3.25) 

3 C 



Let us note that in the LHS of Eqs. ( ^^22] ) and (|]2§) there is a term proportional to 
This might appear as a mismatch of the 0-th order and of the first order in e. This is not 
the case, however, since, in these equations, behaves as an auxiliary field which allows 
to connect these two equations. Indeed, replacing the RHS of Eq. ( 3.22| ) in Eq. ( 3.24| ) one 
recovers the trace of the gravitational field equations of ACDM. 

Finally, from Eq. ( 2.17| ), we find 



^(1,0) / ^ 3^(1,0) / 
= 3^(1,1)'. 



(3.26) 
(3.27) 



At this point, from Eqs. ( p^ ), ( |3l^ ), dOOj ), [substituting the RHS of Eq. 

( 3.22 ) in the RHS of Eq. ( ^.24)) ] and ( 3.26 ), we obtain the perturbation equations in the 
synchronous gauge in general relativity (see Appendix |^ and Refs. [^, |9^). Instead, 
from the other equations we are able to get their correction terms, i.e. 5^^''^\ 2V^^/^*^^'^) + 
(l/3)V2v2x(^'^). In particular, knowing that Ti" = U'Ti + 71^, and after some tedious 
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calculations, we derive the following relations^ 



da 

^(1,0) 



1 

+ 3 



5(1,0) . 



da / 

(3.28) 



3 1 AttoP'p^ 



njl2 njiZ 

ISn^ - 35n' + 23^ - 6^ 



s 

da 



n_ii2 -11/3 \ 



+ 



I— 



(0) 



25-H^ - im' + 6 



/2 



+ 77^-2 



+ 



n 



da 



1^ 

UvraVS^ 

n 



dA 



47ra2p 



(0) 



da / 



(3.29) 



Moreover, by linearizing the solution of the continuity equation ( p. 27 ) we obtain 



V.(^'^)(,y,x) = 4^'^)(x) + \ (<5(i'i)(,y,x) - 4^'^)(x)) . 



(3.30) 



^In order to illustrate better that the background is the same of the ACDM model for <g 1, let us 
consider for example Eq. (3.2^). In the RHS of this equation most of the coefficients are proportional to TH? 



or H' /H^. In this case, they describe the same background as the ACDM model because the corrections 
are of order e^. In other words, they are negligible at first order in e. 
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We denote by a subscript the condition at the present time of the referred quantity. Then 



Iy2y2^(l,l) 

3 



+ isn 



2 [im' - - 17^ + 12— j e 

2 {bh' - m' + 2^) i + WpW (^_,^^ I 5(1,0) 



•T//2 

' 10H' + 2^ 



da 



25^' 



197^' + (-a^] + 3 ( im' 

da . 



12 



7?^ 



I — 

da 



1 



2 JO) 



+ 2-7^ + 2 



n 



-sn^ + m' -—^ 



-a^ 1 + 4( 3-H^ 
da , 



4^' + ^)^ 



127ra2pm 



-im + 2 



S 
da 



+ 4 



+ 



127ra2p; 



(0) 



(1,1) 



(3.31) 



At this point, one can remove the residual gauge ambiguity of the synchronous coordinates 
by imposing that^ '^o^'^^ = ~(l/2)V^Xo^'^^- Therefore, evaluating Eq. ( p. 31 ) to present 
time, we can determine immediately the value of V'o^'^''- 

Eqs. (3.28), ( 3.29| ) and (3.31), which govern the lowest-order modified gravity correc- 
tions w.r.t. ACDM in the behavior of scalar perturbations, represent the main result of this 
paper. Let us stress two important facts: 1) all our expressions are completely determined 
by using the well-known results of the ACDM model (see Appendix ^) and the dynamical 
solution of ^; 2) through this method it is possible to obtain analytically a more precise 
result if we consider the next orders in e. 

In the next section we will analyze in detail the evolution of perturbations in the 
Poisson gauge. 



4. From the synchronous to the Poisson gauge 

In this section we obtain the evolution of cosmological perturbations in the Poisson gauge'' 
(also known as the conformal Newtonian gauge or the longitudinal gauge) through a gauge 
transformation of the results obtained in the synchronous gauge in the previous section. In 
particular we will follow the approach used in Ref. |9^] (e.g. see also [^]). Setting linear 

^In addition, we can conclude that Jq^' = — (l/2)V^Xo^'', Appendix ^ 

^In the Poisson gauge one scalar degree of freedom is eliminated from the goi component of the metric, 
and one scalar and two vector degrees of freedom are eliminated from Qij. 
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vector and tensor modes to zero, the flat linear metric becomes ds^ = a^[— (1 + 2$p^^)(iry^ + 
(1 — 2^^p^)5ijdx^dx^]. Instead, perturbing the mass-density and fluid four-velocity we get 



PSJ{1 + 4 ) and n'^ = {5^ + 4^'')/a, where = -^^^^ and = d' 



,(1) 



Considering the perturbations at the same space-time coordinate values, the syn- 
chronous gauge and the conformal Newtonian gauge can be related by the following rela- 



tions ||93 



2cl>« 

^(1) 



X 



(1) // , 



(1)' 



2^ 



(4.1) 
(4.2) 

(4.3) 

(4.4) 



At this point, as we have already done in section ^, we can split \ '^^\ and 
I'p^'' in the following way 



$(1) 
p 

p 

"p 



,5(i,o) + ,^(i,i) + ._ 

^(1.0)+,^{l,l) + ._ 
5(1.0) +,5(1.1) + 



,(1.0) 



+ ev 



(1,1) 



+ 



(4.5) 
(4.6) 
(4.7) 
(4.8) 



Also in this case the terms «>J,^'°\ 4^'°^ and v'^'^^ are the perturbation terms 

that one obtains in a ACDM model. Now, we want to study and determine the terms to 
order e. Immediately we note that using Eqs. (|1|), and (|3^) we can obtain Wp^'^^ and, 
consequently, 5p''^\ Indeed 



,(1.1) 



e 2 (0) 



n 



12 



/3 



-18^^ + 35^'?^ - 23— + 6 



s 

da 



2 [im"H - b-H' - 17^ + 6^) C - 4^aVS^ (rn - 2^) 



7r\ 
n) 



V-25(1.0) 



25-H^ - im' + 6 



+ 6 (^w - - 2^^ i + ^Tia^p^^i 



V-25(1.0) ' 



+ 



'n!_\ 



5(1,0) +(_«^) 5(1,0) 'I, . (4.9) 



_3 
da 

dA 

""d^J 



where V ^ stands for the inverse of the Laplacian operator; while from Eqs. (4.3) and 
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(I3.28D we get 



P - o (0) 



^(1,0) / : 



Sim 



+ 



(0) 



+ 4 ( + ^ ) C 



Svra^pm 



(0) 



+ 



3^^ 



(0) 

m 



-187^2 ^35?^' -23^ +6 ^4 



da 



/ V da 



S 
da 



2 (l6?^'^ - 5n' - 17^ + 6^) e - 4^aVL°^ " 2^) 



+ 



3H 



(0) 



25^2 _ 19^/ _^ g 



S 
da 



+ 3 ( 10?^' - - 4— ) e + s^aVS^^ 



V-25(1.0) ' . 



da^ 
(4.10) 



Instead from Eqs. (|l0l), (|^) and ( p9| ) we find 



$(1.1) 
p 



12 



3n^ -m' + -a-^ 



dA 



da / 



3^2 _ 5-^' + 2^ ) ^ 



'1^ 



+ 2vra2p(0) 



^d^ 
da 



V-25(1.0) + 



1 



m-2 



n 



de 

I — 
da 



n J 



3 V day 



(4.11) 



Moreover from Eqs. 



^(1.1) 
p 



(U), (|l0|) and (|3^ ) we obtain 



V-25(1.0) + 



1 



7?^ -2 



n 



dA 

''^d^J 



+ 4(-^ + ^)e]v-.(-)'4(-a|).(.o). 



(4.12) 



Finally, another useful quantity is the anisotropic contribution n(i) which is one of 
the parameters that allow to quantify the departure of f{R) gravity from the standard 
ACDM model @. Indeed in our formalism n(i) = ^^p^ - <S>'^^ = e(^'f,^'^^ - ${,^'^^) because 
= $f,^'°) (see appendix 0). Then 



(4.13) 



5. Comparison with the "quasi-static" approximation 



As it is well known, when we consider scales deep inside the Hubble radius, in order to 
derive the equation of matter perturbations approximately, one uses the quasi-static ap- 
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proximation (for details, for example, see ||3^). This section is devoted to recover this ap- 
proximation, as an important example of how our formalism confront with and recover some 
known results in the literature. In particular we want to calculate the Poisson equation 
for these models with the approach studied in this work. Defining = — (^>p^^ + ^'p^^)/2 
from Eqs. and we get 



(1) 

cS 



(5.1) 



Then, knowing that ^^i = Ks + ^^cff and using Eqs. (|^, (|3T7D and (|1|), we 
obtain 



1 

2 L 



(5.2) 



At this point, let us consider in detail Eq. ( |3.23| ). In that scales we can drop the terms with 
temporal derivatives when we compare them with spatial gradients term in Then 
we find 



} {%" - 2W) -2{n' + W) - 



l87raVS^5(i'i)-eVi'0)'. (5.3) 



As we see from Figs. |l| and |2| , [^/(d^/da)] is less of 1/10. Moreover, considering scales 
where the square of wavenumber k is larger than Ti' and TH? , we conclude that — V^^^^'^-* ~ 
[(87ra^/3)pm^'^*-^'"'^^ — ,^'5^^''^^ ']. Consequently, in this case, the additive term proportional 
to becomes negligible with respect to the other terms in Eq. (|5.2|). Then 



(5.4) 



On the other hand, starting from the literature (for example, see Refs. [77, Q) and 
assuming the quasi-static approximation (i.e. when V^|X| ^ "H^jXl and \X'\ < 'H\X\, 
where* X = (^(1), (/.(i) ')> we find V^^^^^ = -47raVm ^^^Vv^- Then we can 

quickly obtain the same result. Indeed 



(5.5) 



However, let us stress that our formalism is designed in order to easily account for effects 
that go beyond the quasi-static approximation, in particular those taking place on scales 
comparable to the horizon size, where time derivatives cannot be neglected. Indeed they 
could be important, for example, for accurate calculation of the Integrated Sachs- Wolfe 
(ISW) effect and the large-scale matter power spectrum (for example, see |103f|). 



*In general, for these models, one also adds the approximation that corrisponds to \ip{d?V / Aip'^)\ 2> 
|(dy/d(p)|. Let us stress that this condition is automatically satisfied through our approach (see section 

Ml). 
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6. Matter Power Spectrum 

In this section we compute the matter power spectrum P{k,a): 

(5(1) (k, a) 5^1) (k', a)) = {27TfsU^ + k')P{k,a) 
where k = |k|. Within our approach, it becomes 

P{k, a) = P(°) {k, a) + eP(i) {k, a) 



(6.1) 



(6.2) 



where P^^\k,a) is the matter power spectrum in the ACDM model. In terms of the 
primordial power spectrum and the transfer function, we have 



P^''\k,a)=Pprim{k)T\k) 



D{a) 

m 



where D is the growing mode of 5^^'^^^ see Appendix Prom Eq. ( [3.28 ) we get 

PW(A;) = 2A{k,a)P^°\k,a) 

with 



(6.3) 



(6.4) 



3^lm{a)A{k, a) 



15 9 
7 - Y^mia) + -r^A(a) 



+ 2 



1 2 



13 

+ 2l-3 + —nm{a) 



-Ma) -2 



2^ 

da 



+ 



/(a)|2 

de 



4-17m(a) + -rJA(a) 



S 
da 



1 - Om(a) + -^A{a) 



da \n 



e(a 

(6.5) 



where / is the growth rate, r2jn(a) is the density parameter of non-relativistic matter 
and ilA(a) is the density parameter of the cosmological constant in a ACDM model, see 
Appendixes |A| and 0. 

Then, defining the initial solutions as in Section ^T^, we can finally derive P{k, a = 1). 
For simplicity, we fix ^(arec) = and (d^/da) (orec) = —5. 10~^ in order to get ^(a = 1) ~ 
-10-2 and (d^/da) (a = 1) ~ -2.43 • 10"^ 

In Fig ^ we show P{k, a = 1) and 



P 



P(0)(fc,a = 1) 



which measures the relative deviation from the matter power spectrum in a ACDM model. 
At this point, let us stress that, knowing 



dl 
dR 



1 + ee 



(6.7) 
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khlMpc-') kh(Mpc"') 

Figure 3: Illustrative plot of P and (AP/P), as a function of k and for a — 1. The lines, from short 
to long dashes, correspond to e = 10^**, 10~^, 10^^ respectively; the black solid line corresponds to 
e = 10-'^. 



and defining Eq. (|3.15|) today, i.e. 



1 ( ^dg 

:i-Oom + J^OA/2) V 'da 



(6.8) 



a=l 



we can compare a generic model in the literature characterized by Bq and Jrq = {df /dR)o 
with our approach (see the discussion in Section |3.1| ) . In particular we make a comparison 
of our results with the Hu-Sawicki model, see for example |^^. In this case, using their 
definitions, we get 



R + f{R) = -2A 



n R"^ 



(6.9) 



where f{R) = f{R) — R and //jo = (df /dR)o in the notation of [^]. Here we have used an 
approximate definition of their Lagrangian, i.e. when this model describes a background 
similar to ACDM. Then e^(a = 1) = fjio and find 



da 



(n + i; 



a=l 



(1 — r?Om + 17oa/2) 
Om — ^^Oa/2) 



+ n 



(6.10) 



Using the initial solutions as in Section |3.1| we see that we are considering a Hu-Sawicki 
model with n = 0.82855. Finally, from Fig |3| we can conclude that our technique work 
for I //JO I < 10~®. Obviously, in order to get a more precise result we have to add e^-order 
contributions. 



7. Conclusions 

In this paper we have introduced a novel method which allows to describe the background 
evolution as well as cosmological perturbations in f{R) modified gravity models, under the 
main assumption that the background evolution is close to ACDM. Here we restricted our 
analysis to linear perturbations, although the method is completely general and can be 
extended to any perturbative order. 
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Let us conclude by adding some comments. As we stressed before, this approach is 
completely different with respect to previous ones. Indeed, through this analytic technique 
we can determine quantitatively the deviation in the behavior of cosmological perturbations 
between a given f{R) model and ACDM reference model. Moreover, our treatment is 
general in that all the results depend only on the initial conditions that characterize the type 
of f{R) model we are studying. Specifically all our expressions are completely determined 
by using the well-known results of the ACDM model (see Appendix ^) and the dynamical 
solution of ^. 

Our method allows to study structure formation in these models from the largest scales, 
of the order of the Hubble horizon, down to scales deeply inside the Hubble radius, without 
employing the so-called "quasi-static" approximation |Q (see also the review and refs. 
therein). This can be very useful as: 1) one can use this approach including all horizon- 
scale corrections to correctly interpret data on very large scales p^. 



100, 101, [103]; 2) at 



the non-linear level, it can be used to get the expressions for the effect of primordial non- 
Gaussianity on the matter density perturbation in an f{R) cosmology, fully accounting for 
the corrections arising on scales comparable with the Hubble radius [104] (e.g. in ACDM 
cosmology see |105|, p6|, WO^, WM, Wm, |TT0t). 



We believe that this technique can be extended to other classes of modified theories 
of gravity to providing information about the evolution of large-scale structures in the 



universe. 
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A. Evolution of first-order perturbations of ACDM model in the syn- 
chronous gauge 

In this section let us consider briefly the evolution, at the linear level, of relativistic per- 
turbations of a ACDM model in the synchronous gauge. Starting from the line-element 
deflned in section ^, in a perfectly homogeneous and isotropic universe, the momentum 
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constraint, the continuity equation and the energy constraint, respectively, give^ 
V'^^^ + = 4^ + = const. , 



2 JO) {(1) _ (1) 



2,11) 



The evolution equation becomes 



X 



(i)" + 2^X^i)' + lvV^ = -2^«. 



(A.l) 
(A.2) 

SvraV^^^^^ • (A-3) 
(A.4) 



An equation only for the scalar mode x^^^ can be obtained by combining together the 
evolution equation and the energy constraint, 

and, consequently, we can get the equation for the linear density fluctuation 

5(i)" + ^<5{i)'_47raVS^<5« = 0. (A.6) 

The equations above have been obtained in whole generality; one could have used in- 
stead the well-known residual gauge ambiguity of the synchronous coordinates (see, e.g., 
Refs.[lll, 112, to simplify their form. For instance, one could fix Xo^^ so that V^Xo^^ = 
—26^\ and thus the x^^^ evolution equation takes the same form as that for 6^^\ i.e. 
6^^^ = — (l/2)V^x^^^- Now, this gauge-condition replaced in Eqs. ( |A.4| ) and ( |A.5| ) yields 

(A.7) 



7^X« ' + W A« + Vx« = -2V.W 



Moreover, from the momentum constraint we get'i/'^^) ='0r-(l/6)V2(xW-x{,^^). Finally 
with such a gauge fixing one obtains 6^^^ oc 'D±, where 'D± represent the the growing (+) 
and decaying (— ) solution of the equation 



In what follows, we shall restrict ourselves to the growing mode. Then 

x(') = 2?+(^?)xS'^ ; and V^'^ = rf^ " \ C^+iv) " 1) V^xi,'^ 
where we have defined 2?+(r/o) = ao = 1. Replacing (A.7) in ( [A.9| ) we find 

3. 



-fir, 



const. 



(A.^ 



(A.9) 



(A.IO) 



'For simplicity, in these appendices, we have substituted the double superscript (1,0) with (1). 
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where = 87raVmV(3^^) andi" f{n^) = dlnP+/dlna. Then 



10 



(1) 



m 



(A.ll) 



It may be convenient to define the gravitational potential today <I>o (see Appendix ^ 
through the relation V2^>o = ^Tialpfls^^^ = -27ra2p[,°^V2xi^^ Then x^^^ = -4$o/(3HgJ7, 
Finally we obtain 



1 + 



2/(0, 



Ora) 



3 r^or 



^0 + 



2 V2$, 



(A.12) 



B. Evolution of first-order perturbations of ACDM model in the Poisson 
gauge 

The goal of this Appendix is to briefly recall the results for the linear perturbations in the 
case of a non-vanishing cosmological A term in the Poisson gauge [116, 117|. 

Starting from the line-element defined in section we note immediately that, at linear 
order, the traceless part of the (i-j)-components of Einstein's equations gives <I>^^^ = ^'^^^ = 
<I>. Its trace gives the evolution equation for the linear scalar potential <I> 



+ m^' + {2%' + 'H^)<^ = . 
Selecting only the growing mode solution one can write 

«>(x,7?) =5(r/)$o(x). 



(B.l) 



(B.2) 



where 4>o is the peculiar gravitational potential, linearly extrapolated to the present time 
rjQ and g{ri) = D^{r])/a{7]) is the so called growth-suppression factor. The exact form of g 
can be found in Refs. |113| , |114 , 115]. In the A = case, g = 1. A very good approximation 
for 5f as a function of redshift z is given in Refs. [113, 114] 



5 oc In^^'' - + (1 + f^m/2) (1 + f^A/70) 



-1 



(B.3) 



with = J7om(i + zf/E\z), = noA/E^z), E{z) = {I + z)n{z)/no = [f^Om(l + 

z)^ + ^Iqa]^^"^ and il-om, ^OA = 1 — ^^Om, the present-day density parameters of non- 
relativistic matter and cosmological constant, respectively. According to our normalization, 
g{z = 0) = 1. The energy and momentum constraints provide the density and velocity 
fluctuations in terms of <I> (see, for example, Ref. |l|] and lll9|, |l20| for the A case) 



,(1) 



(B.4) 
(B.5) 



^'^Specifically /(Sim) = 1 + g' /{T-Lg) where g{z) = "D+ja the growth suppression factor for a ACDM 
Universe (see Appendix H). 
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